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Abstract 

A three-dimensional viscous structured / 
unstructured hybrid scheme has been developed for 
numerical computation of high Reynolds number 
turbomachinery flows. The procedure allows an 
efficient structured solver to be employed in the 
densely clustered, high aspect-ratio grid around the 
viscous regions near solid surfaces, while employing 
an unstructured solver elsewhere in the flow domain 
to add flexibility in mesh generation. Test results for 
an inviscid flow over an external transonic wing and a 
Navier-Stokes flow for an internal annular cascade 
are presented. 

Introduction 

In modern turbomachinery designs, the rotor and 
stator blade rows often possess extreme turning 
angles such that the flow direction deviates greatly 
from the axial direction. These geometries frequently 
cause difficulties in generating structured meshes [1]. 
In order to impose point-to-point periodicity on 
periodic boundaries in the domain, the computational 
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grids are typically forced into highly skewed shapes 
that decrease the accuracy of the solvers. One remedy 
for reducing the grid skewness is to forego the point- 
to-point correspondence on the periodic boundaries 
and to interpolate solution values where necessary. 
The disadvantage of the non-periodic grid approach is 
that the local density of grid lines on either side of a 
periodic boundary may be vastly different. Hence, 
flow features resolved on one side of the domain 
boundary may be diffused or lost when interpolated to 
a coarser grid on the other side. Another solution to 
grid skewness is to use unstructured meshes. 
However, the cost and technology of present day 
unstructured solvers still leave much room for 
improvement. 

The computation of both inviscid and viscous 
flows on unstructured triangular meshes in two 
dimensions, and tetrahedral meshes in three 
dimensions, has matured significantly in recent years 
[2-7]. The ability of unstructured solvers to handle 
complex geometries has proven to be valuable for 
many computations. However, the geometric 
flexibility of unstructured grid solvers is also the 
source of their disadvantages when compared to 
structured solvers. Due to the naturally-ordered data 
connectivity of structured grids, structured solvers 
require much less memory (especially for implicit 
schemes), implicit structured solvers are straight 
forward to code, and turbulence modeling is easier to 
implement. 

The shortcomings of both structured and 
unstructured methodologies are continuously being 
addressed and improved. In the meantime, however, 
it is possible to take advantage of the best properties 
of both methodologies. In the present paper, a 
solution procedure which couples an efficient 
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structured solver with an unstructured solver is 
presented. The hybrid procedure takes advantage of 
the computational efficiency of structured codes and, 
at the same time, is able to benefit from the geometric 
flexibility of unstructured solvers. With the hybrid 
approach, densely packed structured grids can be 
placed in the highly viscous regions near solid 
surfaces, and unstructured grids can be used away 
from solid surfaces [8-10]. This approach can avoid 
the severe grid skewness commonly experienced by 
fully structured grids around turbine blades with high 
turning angles. The procedure has various other 
applications as well, such as providing a means of 
connecting multi-body geometries for Chimera-type 
grids to insure node-to-node correspondence [11-12]. 

The strategy of coupling structured and 
unstructured methods has been implemented for two- 
dimensional turbomachinery computations by many 
researchers in the past [13-15]. The present work 
extends this strategy to three-dimensional 
turbomachinery flows. The hybrid solution technique 
is first tested for an external fixed-wing transonic flow 
case and the result is compared with a structured 
solution for validation. The code is then applied to an 
annular cascade and the result is compared with 
experimental data. 


Formulation 

In this section, the solution procedure for the 
unstructured solver is first described, followed by the 
procedure for the structured solver. Since both the 
structured and the unstructured solvers used in this 
paper have been previously documented, only brief 
outlines of the mathematical and numerical 
formulations are given. Lastly, the coupling procedure 
used to integrate the two solvers is detailed. 


Unstructured Solver 

Governing Equation 

The time dependent, Reynolds averaged, 
compressible Navier-Stokes equations, which express 
the conservation of mass, momentum, and energy, are 
solved. The turbulence viscosity is calculated using 
the high Reynolds number turbulence model of 
Launder and Spalding [16]. The equations of motion, 
written in an integral form for a bounded domain Q 
with a boundary 90, are 


+ JJai F (Q) • AdS (1) 

= Jj^G(Q) • ridS + JJJ q S(QW 

where Q is the unknown vector containing the 
conserved properties 

Q = {p, pu, pv, pw, e Q ,pk, pc} T (2) 

In the above equation, p is the fluid density, u , v y and 
w are the Cartesian velocity components in x, y, and z 
directions, respectively. e 0 is the total energy per unit 
volume. The turbulent kinetic energy and turbulent 
kinetic energy dissipation rate are represented by k 
and & The vector F is the convective flux term, G is 
the viscous term, and S is the source term containing 
the production and destruction of turbulent kinetic 
energy. A complete description of the governing 
equations is presented in Reference [17] and will not 
be detailed here for brevity. 

Inviscid Flux Spatial Discretization 

The inviscid flux across each cell face xris 
computed using Roe’s flux-difference splitting formula 
[18] 

F k = a[f(Q l ) + F(Qfl) - |a|(Q* - Q l )] (3) 

Here, and are the state variables to the left 
and right of the interface k. The matrix A is 
computed by evaluating 


with Roe-averaged quantities so that 

F(Q«)-F(Q l ) = A[Q fl -Qx,] (5) 

is satisfied. 

For a first-order scheme, the primitive variables 
at each cell face are set equal to the cell-centered 
averages on either side of the face. For a higher-order 
scheme, estimation of the state at each face is 
achieved by interpolating the solution at each time 
step with a Taylor series expansion in the 
neighborhood of each cell center. The cell-averaged 
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solution gradients required at the cell center for the 
expansions are computed using Gauss’ theorem by 
evaluating the surface integral for the closed surface 
of the tetrahedra. This process can be simplified using 
geometric invariant features of tetrahedra [4]. The 
resulting second-order formula for the flow state at 
each cell face can be written as 

?A . 3 .3 = 9i + 7 [f (?«, + <3n, + 9« 3 > - 9 n „ ] ( 6 ) 

where the subscripts i%i , n<i , n% denote the nodes 
comprising face f\2 2> °f ce ^ and n 4 corresponds to 
the opposite node. The expansion also requires the 
nodal value of the solution, which can be computed 
from the surrounding cell center data using a second- 
order accurate pseudo-Laplacian averaging procedure 
as suggested by Homes and Connell [19]. The three- 
dimensional extension by Frink [20] is adopted in the 
present calculations. 

The convective terms of the turbulence equations 
are calculated using a first-order accurate scheme in 
the present paper to reduce the computational cost 
and to ensure the numerical stability of the time 
integration [21]. 

Viscous Flux Spatial Discretization 

The evaluation of the viscous term G requires 
first derivatives of the velocity, temperature, and k-e 
values at the cell faces. They are achieved by first 
evaluating the gradient of each required flow quantity 
at the cell center from the known primitive variables 
at each time step. The gradient of the desired 
quantity is obtained by applying the gradient 
theorem, 

(7) 

v Cl 

where Q represents the volume of the domain over 
which the theorem is applied. The scalar quantity 0 
can be the three components of velocity, the 
temperature, or turbulence quantities. In the present 
calculations, the integral domain is defined as the 
individual tetrahedral cell, and the surrounding 
surface area 9G consists of the four triangular 
surfaces covering the cell. This formulation is 
consistent with the numerical procedure of evaluating 
the convective fluxes of the present cell-centered 
scheme. 


Once the gradients of the desired quantities are 
known at the cell center, nodal values are calculated 
using the pseudo-Laplacian averaging mentioned 
earlier for the convective terms. The flux through 
each of the triangular faces in equation (7) is obtained 
by averaging the three nodal values for the triangle. 
Once the gradients of the primitive variables are 
obtained, the shear stresses and can be calculated, 
from which G is evaluated at the cell center. The 
nodal values of these quantities are calculated once 
again by applying the pseudo-Laplacian averaging of 
the surrounding cell center values. The surface flux of 
these quantities in equation (1) is obtained by taking 
the average of the three nodal values for each 
triangular face of each cell. 

Time Integration 

The seven equations of motion are integrated in 
time using an explicit 3-stage Runge-Kutta scheme 
developed by Jameson et al. [22]. The inviscid fluxes 
are evaluated at each time stage using values of 
transport variables obtained at the previous stage of 
the scheme rather then using values from the 
previous iteration [23]. The viscous dissipation and 
the source terms are evaluated prior to the first stage 
and remain constant within that time iteration. 

Local time stepping and implicit residual 
smoothing are used to accelerate the convergence to 
steady state. Inclusion of the viscous terms in the 
residual smoothing procedure was essential to obtain 
a stable and convergent solution. The implicit 
smoothing is also applied to the turbulence equations. 

Structured Solver 

The structured Navier-Stokes analysis in the 
present procedure uses a center-differenced, finite- 
difference scheme. The code was originally a single - 
block external flow solver developed for analysis of 
fixed and rotary wings [24]. Multi-block capability has 
since been added to the code and additional boundary 
condition routines have also been implemented for 
internal flow analysis. 

The structured code solves the unsteady, 
compressible, Reynolds-averaged Navier-Stokes 
equation in strong conservation form for curvilinear 
coordinates 

dQ <?(E-Ev) <?(F - Fv) <?(G - Gv) _ , 0 , 

1 1 1- = u loj 

dx dn dt; 
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Spatial Discretization 


Here, Q is the vector of conserved variables divided by 
the Jacobian of the coordinate transformation. Five 
unknowns are solved instead of seven as in the 
unstructured solver since Baldwin-Lomax algebraic 
[25] turbulence model is used for the structured 
solver. 

Time Integration 


Standard second-order central differencing is 
used for the spatial derivatives. Spectral -radius - 
scaled fourth/second-difference artificial dissipation 
[22] is added for stability and to eliminate oscillations 
near shocks. The viscous terms are evaluated using 
half-point central differencing so that the 
computational stencil for the stress terms uses only 
three nodes in each of the three directions. 


The solution vector is integrated in time using 
the implicit Euler method, which is first-order 
accurate in time. The nonlinear inviscid flux vectors 
are linearized at every time level about their values at 
the previous time level using Taylor expansions, e.g. 

E n+1 = E n + — " AQ n+1 + ofAr 2 ) 

dQ V ) ( 9 ) 

= E” +A"AQ' l+1 


The viscous terms are evaluated explicitly. Explicit 
treatment of the viscous terms still permits the use of 
large time steps since the Reynolds numbers of 
interest here are fairly large. To further reduce 
computational time and memory, the radial/spanwise 
flux derivatives are treated explicitly at the old time 
level[26], but the new values are incorporated as soon 
as they become available. This explicit treatment of 
the span wise flux terms enables the scheme to solve 
the three-dimensional equations by solving one 
spanwise station implicitly at a time. To eliminate 
any dependency the solution may have on the 
spanwise marching direction, the solver reverses the 
marching direction with every spanwise sweep. 

The resulting left-hand-side of the matrix 
equation is approximately factored into alternating 
directions 


(i + At 5|A n )(l + At <5 ? C n )AQ 


/l+l 


= -At 


(E“ -Ev n ' n+1 ), +(F-Fv)J 
+ (G” - Gv n,n+1 ) 


71+1 


( 10 ) 


The implicit operating matrices are then diagonalized 
using a similarity transformation [27]. 


Hybrid Coupling 

Since the structured code has multi-block 
capability, the coupling procedure treats the 
unstructured portion of the hybrid meshes as a 
separate block. The blocks are solved independently 
with respect to each other. Explicit boundary 
conditions are updated at the end of each time 
iteration. The two codes are loosely coupled, 
essentially running independently of each other. The 
only interaction between the structured solver and 
the unstructured solver is through boundary 
conditions. The advantage of this formulation is that 
almost any two given codes can be coupled together. 

In this paper, a central -differencing, finite-difference 
procedure and an upwind cell-centered finite volume 
procedure are coupled. The structured solver is 
formulated in an inertial frame of reference where as 
the unstructured solver is in a blade-fixed coordinate 
system. For non-moving-boundary cases, such as 
pressure-driven turbine stator flows, there is no 
difference between the two formulations. However, for 
moving-boundary cases, care must be taken when 
transferring data between the two solvers depending 
on how boundary conditions are prescribed. For 
example, consider a wing traveling at a constant 
velocity. Since the structured solver assumes the wing 
is moving and the unstructured solver assumes the 
domain is fixed in space, the total energy must be 
converted between the two frames of reference. At the 
end of each iteration, mass, momentum, and energy 
values are communicated between the structured and 
unstructured solvers. 

Since the unstructured solver is a cell-centered 
finite-volume scheme and the structured solver is a 
vertex-based finite-difference scheme, the 
unstructured and the structured grids are overlap by 
one layer at the boundary in order to minimize 
interpolation and extrapolation errors at the 
boundary. By overlapping the grids instead of 
abutting them against each other, the boundary nodes 
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for the structured solver are interpolated from the 
interior cell center values of the unstructured solver. 
Similarly, the cell boundary face values of the 
unstructured solver are interpolated from interior 
structured nodes. An illustration of the overlapped 
grid in 2D is shown in Fig. 1. 



Figure 1. Structured-unstructured mesh overlapping 
by one layer. 

Once the grid construction is completed, a table is 
generated to index unstructured nodes to their 
corresponding structured nodes, and vise-versa. The 
index is similar to the table generated between 
structured block boundaries for multi-block 
calculations. 

Results and Discussion 

Validation for both the structured [24] and the 
unstructured solvers [28, 17] have been documented 
previously and will not be presented here. To test the 
structured-unstructured hybrid analysis, calculations 
for a transonic inviscid flow over a fixed wing and a 
viscous flow for a turbine stator blade are presented. 

For all calculations in this paper, a structured 
grid is first generated for the entire geometry. The 
unstructured grids are generated by subdividing a 
portion of the structured cells each into six 
unstructured tetrahedra (Fig. 2). This procedure 
simplifies grid generation and is suitable for the 
purpose of code validation. 



Figure 2. Unstructured tetrahedral cells generated 
using structured cells. 

Transonic Wing 

The first calculation selected for the hybrid 
scheme is a transonic wing at a freestream mach 
number of 0.8 and at 1.25° angle-of-attack. For this 
test case, flow periodicity is imposed on the two end- 
wall boundaries to simulate a two-dimensional flow. 

A purely structured solution for this case is first 
obtained for comparison. An O-H grid topology is 
used; that is, an O-grid is generated at each spanwise 
station. For the structured calculation, each spanwise 
station consists of 121 (streamwise) x 41 (normal) 
points. Six spanwise stations are required due to the 
memory allocation algorithms for the multi-block 
structured solver. The number of spanwise stations 
used does not alter the results since the code 
simulates a two-dimensional problem. 

For the structured-unstructured coupled 
procedure, the grid is divided into three zones. The 
inner zone, which wraps around the solid surface 
(121x12), is calculated using the structured solver. 
The middle zone (121x6), which wraps around the 
inner structured zone, is subdivided into tetrahedra 
for the unstructured solver. The outer zone, which 
extends from the middle zone to the free stream, is 
once again calculated using the structured solver. 

Fig. 3 illustrates the three-zones described for this 
case. 
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Figure 3. Structured-unstructured zones for the 
external wing calculation. 


The actual grid used for the calculation is shown in 
Fig. 4. 



Figure 4. Structured-unstructured grid for the 2-D 
wing calculation. 

The purpose of dividing the calculation into three 
zones is to simulate a Chimera type of grid where a 
simple background grid is generated for the entire 
domain and a structured grid is used around a solid 
object. Instead of interpolating between the two grids, 
an unstructured solver can be placed to connect the 
two set of grids and obtain node-to-node 
correspondence which remove the needs for 
interpolation. For a finite volume approach, this also 
insures flow conservation. This approach is also 
recommended for multi-body calculations. Many 
researchers have demonstrated structured- 


unstructured coupling for multi-element airfoil 
calculations where structured grids are place near the 
solid surfaces and an unstructured solver connects the 
structured grids and extends to the farfield. Since 
structured solvers are generally more efficient, it 
makes sense to minimize unstructured zones. 

The unstructured solver can be run in both first- 
and second-order modes, while the structured solver 
is always second-order accurate. Surface pressure 
distributions for both first- and second-order hybrid 
solutions, as well as the fully structured solution are 
shown in Fig 5. Both hybrid solutions agree very 
closely with the standard structured solution, but the 
first-order unstructured solution exhibits a larger 
deviation on the lower surface. 



Figure 5. Surface pressure distribution. 

To examine the coupling between the structured- 
unstructured solver across the three-zones, pressure 
contour lines are shown in Fig. 6 for the first-order 
and the second-order coupling, along with the single 
block structured solver for comparison. For the first- 
order coupling, a slight discontinuity exists across the 
zonal boundary. However, the discontinuity is small 
and barely perceptible on a global scale. For the 
second-order coupling, the contours show the flow is 
smooth with little indication of the presence of the 
zonal boundaries. The mass flows in and out of the 
sandwiched unstructured zone are within 0.03% of 
each other, indicating that good flow conservation is 
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achieved. The unstructured shock is somewhat more 
compact then the structured solver. One reason for 
this is that every structured cell becomes six 
tetrahedra so that in the stream wise direction, the 
unstructured cells are twice as dense as the 
structured cells. The present hybrid procedure 
converges considerably slower then the single 
structured grid due to the explicit formulation of the 
incorporated unstructured solver, which requires a 
low CFL number. 

3-D Viscous Flow Through Turbine Blades- 
Hybrid Procedure 

The second test case selected is a three- 
dimensional annular cascade [29]. The geometry 
consists of an annular ring of 36 turbine stator vanes. 
The blades are 38.10 mm in span, untwisted, and of 
constant profile with an axial chord of 38.23 mm. The 
stator has a tip diameter of 508 mm and a hub-to-tip 
radius ratio of 0.85. The inlet flow angle is parallel to 
the axis of the cascade. The Reynolds number based 
on the inlet total quantities and axial chord length is 
898,650. The inlet total pressure and total 
temperature at one axial-chord length upstream of 
the blade leading edge are known from experiment. 
The exit hub static pressure to inlet total pressure 
ratio is 0.65 at 2.6 axial-chord lengths down stream of 
the blade trailing edge. 

For the hybrid calculation, the computational 
domain extends one axial-chord length upstream of 
the blade leading edge and 2.6 axial-chord lengths 
down stream of the blade trailing edge. The exit hub 
static pressure to inlet total pressure ratio is known 
from experiment. The grid has 105 (inlet to exit) x29 
(blade to blade) xl5 (hub to tip) nodes and 23,000 cells 
are subdivided into unstructured cells. The structured 
grid wraps around the solid surfaces of the hub and 
the turbine blade while the rest of the domain is filled 
with unstructured cells. The hybrid grid at one radial 
station for the annular cascade is shown in Fig. 7. For 
the present grid, only 15 hub to shroud stations are 
used and the grid density near the trailing edge is 
very coarse (only 3 points defining the rounded 
trailing edge). While such a coarse grid is not 
acceptable for accurate prediction of viscous effects, it 
is sufficient for monitoring the communication 
between the structured and unstructured zones and to 
predict the overall flow field. 

A no-slip boundary condition is applied at the hub 
and on the blade surfaces, and an inviscid slip 
condition is applied on the shroud. In the coupled 
calculation, the unstructured solver is run without 


turbulence modeling because the effects of turbulence 
are much reduced outside the boundary layer. 

The pressure contour lines across the structured 
and unstructured zones are plotted in Fig. 8. The 
same degree of smoothness as in the previous wing 
case is observed. In Fig. 9, the chord wise surface 
static pressures normalized by the inlet total pressure 
are compared with experiment at 13.3%, 50%, and 
86.7% span. An overall agreement is found accept 
near the trailing edge where the computation on the 
present coarse grid did not pick up the peak values. 


Concluding Remarks 

A three-dimensional, unstructured, Navier- 
Stokes flow solver has been coupled with a three- 
dimensional structured code to allow structured- 
unstructured hybrid calculations. The two codes are 
loosely coupled and interact only through boundary 
conditions. The structured solver is a center- 
differenced finite-difference scheme and the 
unstructured solver is an upwind-differenced finite- 
volume scheme. 

The hybrid procedure has been tested for a 
transonic wing and an annular cascade, and good 
results have been obtained. The flow properties across 
the boundary between the structured and 
unstructured portions of the grid are smooth and well 
behaved. The results show that two distinct 
discretization techniques can be coupled together with 
little effect on the solution accuracy, as long as both 
are of same order. 
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first-order coupling 


second-order coupling 





Figure 6. Pressure contour lines across the structured-unstructured zonal interface 





Figure 7. Structured-unstructured hybrid grid used Figure 8. Pressure contour lines across the 
for the annular cascade flow calculation. structured-unstructured zonal interface. 
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Figure 9. Static pressure distribution on the blade at 13.3% (a), 50.0% (b), and 86.7% (c) span. 
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